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("The  Defect  Correction  Principle  and  Discretization  Methods")  published 
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appear  soon  In  Numerlsche  Mathematlk.  In  this  report  many  useful  Insights  are 
furnished  Into  the  historical  development  of  estimation  of  error  of  discretiza- 
tion methods.  In  particular,  Stetter's  work  appears  to  constitute  a significant 
advance  over  Llndberg's  work,  referred  to  In  this  report,  on  asymptotic  error 
estimates. 
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SECTION  I 
INTRODUCTION 

There  Is  a serious  uncertainty  about  the  true  accuracy  of  certain  numerical 
schemes  employed  for  obtaining  approximate  solutions  of  systems  of  nonlinear 
PDEs  (partial  differential  equations),  such  as  In  the  hydrocodes  and  wavecodes. 
Until  close  error  estimates  are  proved,  the  output  of  any  numerical  scheme  Is 
open  to  the  criticism  of  being  little  more  than  a guess  about  the  exact  solution. 
Furthermore,  even  the  existence  of  solutions  Is  often  In  question. 

In  this  report,  two  error  estimation  procedures  are  Introduced  and 
Illustrated.  One  procedure  Is  called  the  asymptotic  error  estimation  procedure 
(ref.  1),  and  the  other  Is  called  the  error  gradient  procedure  (ref.  2).  Also 
Included  In  this  report  are  the  results  of  some  computations  which  were  carried 
out  with  the  Intention  of  providing  concrete  Illustrations  of  the  two  procedures. 
Four  ODEs  (ordinary  differential  equations)  are  considered  In  some  detail.  It 
Is  Intended  to  present  the  results  of  corresponding  computations  for  PDEs  In  a 
later  report. 

It  Is  hoped  that  this  report  will  render  a valuable  Instructional  function, 
namely,  to  assist  numerical  analysts  In  obtaining  accurate  estimates  of  the 
errors  Incurred  by  various  computational  schemes,  so  that  an  optimum  decision 
among  the  various  available  options  may  be  made. 

The  authors  also  hope  that  the  specific  computations  which  are  presented 
herein  will  serve  to  stimulate  further  work  along  these  lines.  They  will 
welcome  any  comments  as  well  as  Information  concerning  further  computational 
work  which  has  been  undertaken. 


1.  Llndberg,  B.,  Error  Estimation  and  Iterative  Improvements  for  the  Numerical 
Solution  of  Operator  Equations.  Report  No.  VIVCDCS-R-76-82fl.  nppartmpnt  nf 
Computer  Sciences,  University  of  Illinois,  Champaign,  Illinois,  July  1976. 

2.  Hicks,  D.,  A Procedure  to  Produce  Estimates  for  Difference  Approximations 
to  Differential  Equations,  I.  Introduction  and  Illustration  of  the  Basic 
Concejjt,  Report  No.  SAND-75-D487,  Sandla  Laboratories,  Klrtland  Air  Force 
Base,  New  Mexico,  January  1976. 
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SECTION  II 

THE  ASYMPTOTIC  ERROR  ESTIMATION  PROCEDURE 

1 . ILLUSTRATION  ON  ALGEBRAIC  EQUATIONS  AND  ODEs 

In  this  section  we  Illustrate  the  asymptotic  error  estimation  procedure  on 
two  examples.  First,  to  bring  out  the  basic  Ideas  as  clearly  as  possible,  the 
comparatively  simple  problem  of  locating  the  zero  of  a real -valued  function  of 
a real  variable  Is  considered.  Than  an  ODE  problem  Is  discussed.  In  the 
following  subsection  a PDE  problem  Is  discussed  briefly,  and  then  the  theoreti- 
cal foundations  are  presented. 

a.  Exampl e #1 . 

Let  F:  R -*■  R;  that  Is,  F Is  a real -valued  function  defined  on  the 
whole  real  number  system.  For  each  sufficiently  small  positive  number  h (hence- 
forth, all  h means  all  h In  some  fixed  Interval  0 < h < hQ)  let  the  family  of 
functions  <j>h : R + R be  a p-th  order  approximation  to  F In  the  sense  that  for 
all  real  numbers  x and  all  h the  equality  (x)  * F(x)  + 0(hp)  holds.  As 
usual,  0(hp)  denotes  a quantity  depending  on  x and  h whose  absolute  value  Is 
<_  C hp,  where  C Is  a positive  number  Independent  of  x and  h,  and  p Is  a fixed 
positive  Integer.  (Here  the  parameter  h Is  somewhat  artificial,  but  In 
practical  applications  It  will  play  the  role  of  a mesh- Increment  In  a finite- 
difference  scheme.)  We  make  the  additional  hypothesis  that  F and  are  twice 
continuously  and  boundedly  differentiable;  l.e.,  there  exists  a positive  con- 
stant A such  that  | F* (x ) | , |F"(x)|,  Uh'(x)|,  and  |$h"(x)|  <_  A for  all  x and  h. 
Also,  we  assume  that  ♦h'(x)  ■ F'(x)  + 0(hp),  where  0(hp)  Is  as  defined  above. 

(Of  course,  two  or  more  appearances  of  the  expression  0(hp)  refer  In  general 
to  quite  different  functions  of  x and  h.) 

As  Indicated  previously,  we  are  interested  here  In  Illustrating  how  the 
asymptotic  error  estimation  procedure  may  be  used  to  estimate  a zero  of  the 
function  F.  Suppose  that  the  equations  F(x)  * 0 and  th(x)  ■ 0 both  possess 
unique  solutions,  denoted  y and  n(h)  respectively,  and  that  It  Is  easier  to 
determine  n(h)  (for  all  h)  than  y;  what  can  be  said  about  the  accuracy  with 
which  n(h)  approximates  y? 


The  quantity  $h(y)  is  termed  the  local  discretization  error;  note  that 
♦h(y)  ■ *h(y)  - F(y),  which  by  the  original  hypothesis  Is  0(hp).  However, 
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this  does  not  Imply  that  n(h)  Is  close  to  y;  roughly  speaking,  such  an  undesir- 
able situation  may  occur  when  the  graphs  of  F and  $h  are  quite  flat  In  the 
vicinity  of  y.  For  example,  If  F(x)  3 x3  and  <l>h(x)  * x3  - h (so  that  p * 1), 

then  y * 0 and  n(h)  3 h1/3,  so  that  n(h)  - y Is  large  In  comparison  with  h. 

Note  that.  In  this  case,  F'(y)  3 0,  In  order  to  make  the  general  situation 
quite  clear,  we  Introduce  the  following  concept  of  stability. 

Let  z be  any  real  number;  since  |<f»h'(x)|  <_  A for  all  x and  h, 

Uh(£)  - <frh(z)I  < r whenever  |£  - z|  < r/A,  r being  any  specified  positive 

number.  This  shows  that,  given  r,  the  set 

>c£(z;r)  3 js:  Uh(c)  - 4>h( z ) | < r for  all  h| 

contains  the  Interval  |£  - z|  < r/A.  We  now  define 

R 3 R (z;r)  3 min  |4(|,U)|(  £ In  (z,r),  0 < h < hQ 

We  now  say  that  the  family  of  functions  Is  stable  at  z provided  that  there 
exist  positive  constants  S and  r such  that,  for  all  h,  the  Inequalities 

I^U1)  - *h(z)|  < r,  |*h(£2)  - *h(z)|  < r 

Imply  that 

1C1  - £2|  < S |#h(5l)  - ♦h(€2) | 

Now,  If  for  some  choice  of  r we  find  that  R,  as  defined  above.  Is  strictly 
positive,  then  by  choosing  S 3 R"1  we  see  that  the  family  <j>h  Is  stable  at  z. 

Roughly  stated,  the  family  of  functions  Is  stable  at  z If  these 
functions  do  not  have  a flat  spot  at  z.  If  these  functions  have  a flat  spot 
at  y then  the  numbers  n(h)  may  differ  significantly  from  y. 

Now  suppose  that  for  a certain  function  z(h)  there  exist  positive 
Integers  p and  M,  with  p £ M,  such  that 

M 

z(h)  * ^ h^e..  Independent  of  h 
J-P 


7 


n(h)  * y + V hJei  + °(hM+1) 
J-P 


then  the  right  side  of  this  last  equation  Is  termed  a p-M  asymptotic  expansion 
of  n(h) . 

Now,  suppose  that  a p-M  approximation,  z(h),  to  y Is  available,  with 
M+1  <_  2p.  (Recall  that  p <^M.)  Then  It  Is  easily  shown  that 

♦h(y)  + F(n)  - o(hM+1)  (1) 

This  Is  proven  as  follows:  By  Taylor's  theorem 

*h(n)  3 *h[y  + 2 + °(hM+1)j  3 ♦h(y)  + *h(y)|*  + °(hM+1)j  + " [z  + °(hM+1)J 

where  y Is  some  suitably  chosen  number  between  y and  n.  Since  |$h"(y)|  <_  A, 
♦j,'(y)  * F'(y)  + 0(hp),  and  2p  >_M+1,  It  Is  now  evident  that  the  above  equation 
simplifies  to 

♦h(n)  » <Ph(y)  + F'(y)z  + o(hM+1) 

or,  since  $h(n)  a 0, 

♦h(y)  - - F'(y)z  + o(hM+1)  (2) 

An  entirely  similar  argument  furnishes  the  equality  F(n)  * F(y)  + 

F'(y)z  + 0(hM+1),  and  since  F(y)  - 0,  we  obtain 

F(n)  - F'(y)z  + o(hM+1)  (3) 

Addition  of  equation  2 and  equation  3 furnishes  the  desired  equality  of 
equation  1 . 

Thus,  the  local  discretization  error  of  the  approximation  + F(n) 

Is  of  higher  order  than  the  local  discretization  error  of  $.(•)»  and  this  fact 

E n 

suggests  that  the  solution  n of  the  equation 
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$h(nE)  + F(n)  = 0 

Is  closer  to  y than  n is.  This  Is  in  fact  true  (asymptotically)  under  the 
additional  hypothesis  that  <|>h  is  stable.  This  is  demonstrated  as  follows: 

♦ h^)  - ♦ft(y)  ■ " [+h(y)  + F(n)]  * o(hM+1) 

Therefore,  for  sufficiently  small  h,  |<f>h(nE)  - <j>h(y)|  < r,  and  now  by  the 
stability  hypothesis  we  obtain 

|nE  - y|  <_  S Uh(nE)  - 4>h(y ) | = o(hM+1) 


nE  = y + o(hM+1) 

which  shows  that  nE  constitutes  an  Improvement  over  n. 

We  can  use  this  last  result  to  estimate  asymptotically  the  error  n - y. 


n - y = r,  - nE  + o(hM+1) 

This  simple  observation  Is  the  basis  of  the  asymptotic  error  estimation 
procedure. 


Observe  that  if 


<{>E  * F + o(hM+1) 


♦h(y>  + ♦h 3 °(hM+1) 


and  therefore  the  solution  rf^  to 


♦h(^)  + ♦h^T) 


9 


Let  E be  the  Banach  space  C1  [o,  1 3 » the  family  of  all  real-valued 
functions  possessing  a continuous  first  derivative  on  the  closed  Interval 
[0,  l],  with  norm  defined  by 

I |y| lE  3 |y(0) I + max  | y ' (x) | 

0 < x < 1 

The  Banach  space  E°  will  be  chosen  as  R x C [0,  l];  that  Is,  the  set  of 
all  ordered  pairs  where  a is  a real  number  and  f is  a continuous  real-valued 
function  defined  on  [0,  l],  the  norm  being  defined  by 

I Ifi  lEo  3 l<*l  + max  | f (x) | 

0 < x < 1 


The  operator  F:E  -*■  E°  will  be  defined  by 


Then,  clearly,  the  aforementioned  initial-value  problem  can  be  restated  in  the 
form  F(y)  = 0. 

We  discretize  the  problem  In  the  following  manner:  For  any  positive 
Integer  N let  h * 1/N,  let  Eh  be  the  Banach  space  consisting  of  all  (N+l)-tuples 
(ordered  sets  of  N+l  numbers)  with  norm 

y _ y I 

1 1 y0»  yN  1 1 3 |y0|  + max  -■^*1h — 

o < j < N-l 


10 
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and  let  E be  projected  onto  E^  by  the  obvious  mapping 

Ah(y)  * [y(0),  y(h)-  y(2h) , • • -J 
Similarly,  let  E^  consist  of  all  objects  of  the  form 


with  norm  |a|  + max  |Cj|,  while  E°  Is  projected  onto  E^  by  the  mapping 


Next,  the  operator  F:E  -*•  E°  Is  replaced  by  the  operator  4>h:E°  -*■  E°  as  follows: 
For  y « (yQ,  yx  , y2,*«*). 


Finally,  the  higher-order  discrete  approximation  4^  Is  defined  by 

’ \ 

y 2 - yi  yi  + y2  I 

— fi 2 — •••*/ 

Note  that  Is  accurate  to  second  order,  while  4>h  Is  only  first-order  accurate. 

To  sum  up,  we  list  step-by-step  the  procedure  to  be  employed  In  deriving 
an  asymptotic  error  estimate,  and  then  Illustrate  the  procedure  In  the  case  of 
■the  Initial-value  problem  described  In  Example  #2. 

Step  1 : Construct  a discrete  approximation  to  the  operator  F. 

Step  2:  Construct  a second  approximation  ^ to  F such  that  4.^  Is  of 
higher  order  accuracy  than  4>h. 

Step  3:  Find  n such  that  4>h(n)  • 0. 

Step  4:  Find  nE  such  that  4>h(nE)  ■ - 4>E(n). 

Step  5:  The  error  estimate  Is  n - nE. 


( 


y\  - y0  *o  + 
i 2 — 


1 
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Referring  to  Example  #2,  the  steps  assume  the  following  form: 
Step  1 : 


(Note  that  <j>h  Is  of  the  first  order.) 

Step  2: 

= ^ - ^1+12  0 < j < N - 1^ 

(<{>E  Is  of  second  order.) 


Step  3: 

The  solution  to 


Is 

Step  4: 

The  solution  to 

Is 


(n)  = 0 

rij  - (1  + h)j 

♦h  (^E)  + ♦h  (n)  * 0 

-5  ■ «’  * h>J  (’  * &*) 


Step  5: 

The  error  estimate  Is 


(n  - nE)  J • - ^ (1  + h)J_1 

while  the  actual  error  Is  (1  + h)^  - e^. 

Although  the  next  report  contains  a considerable  amount  of  numerical 
results.  It  appears  worthwhile  to  present  here  a few  numbers  Illustrating  the 
problem  under  consideration.  If  N ■ 20,  h ■ 1/20,  and  j ■ 20,  then 


J 
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(n  - nE)  j * - (1.05)19  » - 0.063  (error  estimate) 


(1  + h)^  - e^  = (1  .05)20  - e w - 0.065  (actual  error) 


Note  that  nj  « 2.716,  so  that  the  Improved  error  estimate,  nj 
- 0.002,  In  contrast  to  the  much  cruder  value  of  - 0.065  obtained  with 

2.  ILLUSTRATION  ON  PDEs 


Is 


Consider  the  simple  problem 


ff  + §f  « 0.  w(0,x)  - f(x) 

The  solution  Is  given,  of  course,  by  w(t,x)  a f(x  - t). 

Now,  for  each  function  w(t,x)  we  define  the  operator  F by 


F(w) 


w(0,x)  - f(x) 

3W  3W 

3t  3x 


(We  omit  the  detailed  definition  of  the  spaces  E and  E°  and  of  their  norms.) 

Let  and  <|>jj  be  discrete  approximations  to  F with  <j>j^  of  higher  order  accuracy 
than  <i>h.  Let  v be  the  solution  of  4>h(v)  * 0 and  v^  be  the  solution  to 
$h  (v£)  + (v)  = 0.  Then  v - v£  Is  the  asymptotic  error  estimate,  that  Is, 

the  approximation  to  v - w furnished  by  the  asymptotic  error  estimation  pro- 
cedure. 

For  example,  let  <|>h  be  determined  by  taking  as  the  difference  approximation 
to  the  POE  the  equation 


that  Is, 


vn+1  - vn  vn  - vn 

vj + j V 


At 


AX 


0 


♦h(v) 


n+i 


At 
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Note  that  *h  Is  accurate  of  order  0(At  + ax).  Next,  let  «j>^  be  determined  by  a 
more  accurate  scheme,  say  the  LWR  method  (Rlchtmyer's  version  of  the  Lax- 
Wendroff  scheme).  That  Is  (ref.  3,  p.  303), 


Note  that  this  scheme  Is  of  order  0(at2  + ax2). 

3.  THEORY  OF  ASYMPTOTIC  ERROR  ESTIMATION 

Definition:  Is  said  to  be  a p-M  approximation  to  the  operator  F provided 

there  exists  an  operator  F on  E which  Is  of  order  hp  and  satisfies,  for  all 
vectors  2 In  E, 

♦h  (4h  2)'  {F(2>  + f,<2>  } * o(hM*‘) 

Furthermore,  If  for  each  vector  z there  exist  vectors  fp,**»,  fM  In  E°  (Indepen- 
dent of  h)  such  that 

f(z)  • 2 "V  fV 
v*p 

then  we  say  that  <frh  - F admits  an  asymptotic  expansion  to  the  order  M.  Finally, 
$h  Is  said  to  be  a Llpschltzlan  p-M  approximation  to  F provided  there  exist 
positive  constants  L and  b such  that,  for  all  h, 

\\F(yX)  - F(y2)  ||  < hp  L ||yl  - y2|| 

whenever  My1  - y| | < b and  ||y2  - y| | < b.  (Recall  that  y Is  the  unique 
solution  to  F(y)  * 0.) 


3.  Rlchtmyer,  R.  D.,  and  Morton,  K.  W.,  Difference  Methods  for  Initial  Value 
Problems.  Interscience,  New  York,  1967!  
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We  Illustrate  the  above  definition  on  Example  #2.  Note  that 


* hL"  - z'(x)  + z"(x)5r  + z"'(x)£  + 0(h3) 


Therefore,  In  this  case. 


F(z)  • ^ hv  fv  (z)  - z"(x)yt-  + z"'(x){£ 

V*1 

where 

fj(z)  = (f~)  and  f2(z)  * (§— ) 

Now  let  us  consider  under  what  conditions  <f>h  Is  a Llpschltzlan  1-1  approximation 
to  F In  Example  #2.  Recall  that  y = ex  and  that 


I ly1  - y2l I 3 ly1 (0)  - y 2 ( 0 ) I + max  I Qy1 (x ) - y2(x)]' I 


0<  x <1 


Thus,  for  M a 1 and  F - fj. 


l|F(yl)  - F(y2)  ||  * M f1(y2)  - f^y2)  ||  * max 


.IV  . fv2Vl 


Thus,  If  s > 2,  then  <ph  Is  a Llpschltzlan  1-1  approximation  to  F.  Furthermore, 
if  s >_  3,  then  ^ Is  a Llpschltzlan  1-2  approximation  to  F. 

Definition : ^ Is  said  to  be  h*^  - Llpschltzlan  at  £ provided  there  exist 
positive  constants  L and  d such  that  the  pair  of  inequalities 

lie2  - £ 1 1 id,  Ms2  - til  < d 


Imply  that 


♦hU1)  - *h(£2)  1 1 < L 1 1 t1  - C2  1 1 h* 
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I 


Therefore,  In  this  case,  <|>h  Is  h"1  - Llpschltzlan  at  every  e. 

Definition;  The  global  discretization  error  Is  defined  as 

e(h)  ■ n(h)  - a^y 

In  Example  #2,  n(h)i  * (1  + h)f  and  a^  • ex1;  therefore, 

e(h)1  ■ (1  ♦ h)1  - ex1 

Definition:  n(h)  Is  said  to  be  convergent  to  y of  ;>rder  p >_  1 provided 
that  e ( h ) * 0(hp). 

In  Example  #2,  one  can  readily  show  that 


|e(h),j | < hxi  ex1 


and  so,  In  this  case,  p • 1. 


Definition;  We  say  that  n Is  a p-M  approximation  to  y provided  there  exist 
M > p and  a vector  z - 0(hp)  such  that  n » y ♦ z ♦ o(hM41).  Moreover,  we  say 
that  the  global  discretization  error  e * n - y admits  an  asymptotic  expansion 
to  the  order  M If  there  exist  vectors  ep,  ep+1,"',  eM  (Independent  of  h)  in  E 
such  that 

M 

z - V'  hve 


V*P 


16 


■a 


— 

j 

AFWL-TR-78-119 

In  Example  #2  one  can  show  that 

*1 

x e h 

e(h)1  - - J-T-  + 0(h2) 

Therefore,  in  this  case,  ex  * - xex/2. 

We  now  proceed  to  state  three  lemmas,  together  with  Indications  of  their 
proofs,  which  will  be  helpful  In  formulating  Theorem  1, 

Lemma  1 . If 

(1)  <t>h  Is  stable  at  c,  and 

(2)  *hU)  + E(n)  - o(hM+1)  (with  M > 0) 

then  the  solution  nE  of  ^(n)  * 0 satisfies 

nE  = C + o(hM+1) 

Indication  of  Proof;  Note  that 

1 1 ^ h (n E ) " Vc)  11*11  *hU)  + *E(n)  || 

By  hypothesis  (2),  1 1 4>h  ( nE  ) - $hU)  1 1 < r whenever  h Is  sufficiently  small. 

Now  by  hypothesis  (1)  we  have 

llnE  - cl  | < s 1 1 $h(nE ) - <t>h(c)||  * s 1 1 ♦h(e)  + <^E(n)|  | * o(hM+1) 

which  Is  equivalent  to  the  assertion  of  the  lemma. 

We  now  Illustrate  this  lemma  In  the  case  of  Example  #2.  Let  s * A^y  (where, 
as  always,  F(y)  « 0)  and  let  4»h(n)  * 0.  Then  Uhy)  + ♦Jj(n)  * 0(h2),  so  that 
M » 1.  Now  A^y^  * ex1,  * (1  + h)*,  so  that  n^,  the  solution  to 
$h(n^)  + $jj(n)  » 0 Is  given  by 

-J-O*")1 

We  therefore  obtain 

nE  ■ ex1  + 0(h2) 


17 
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In  contrast  to  the  weaker  result 

ni  * ex1  + 0(h) 

Remark:  The  purpose  of  Lemma  2 Is  to  provide  conditions  under  which 
hypothesis  (2)  of  Lemma  1 holds. 

Lemma  2.  Assume  that 

(1)  y and  n are  the  unique  solutions  to  F(y)  * 0 and  <|>h(n)  = 0 where  n Is 
a p-M+1  approximation  to  y with  1 < p < M. 

(2)  Is  a consistent  approximation  to  F of  order  M+l  with  p ^ M+l  <_  2p. 

(3)  <|>h  and  <j>^  are  h"1  - Llpschltzlan  at  y. 

(4)  <f>h  Is  a Llpschltzlan  p-M  approximation  to  F. 

Then, 

*h(v) = • *h(n)  + °(hM+1) 

Indication  of  Proof:  In  relation  to  hypothesis  (1)  we  introduce  the 

notation 

n a Ah  (y  + z)  + 6 

where 

z * and  5 * o(hM+2) 

In  relation  to  assumption  (4)  we  Introduce  the  notation 

♦h(v)  • < * ^<y)]  ♦ o (hM*‘) 

where 

F{y)  • o(hp) 


- - - 
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From  assumption  (1)  we  have 

♦h(v)  + ♦h<">  ■ *h(v)  - »h<n>  * ♦h<"> 

■ »h(v)  • ♦h[V»  ♦ z>  ♦ «]  + *h<n>  [Ah(y 

From  assumptions  (1)  and  (3)  the  last  line  may  be  rewritten 

’ *h(v)  - +h(fh(y  * z>]*  *h[ahty  * z>]+  °(h"*1) 

From  assumptions  (2)  and  (4),  furthermore,  this  becomes 

3 aJ  | [F(y)  + F(y)]  - [F(y  + z)  + F(y  + z)]  + F(y  + z)j 
+ 0(hM+1) 

From  assumption  (1),  we  obtain 

3 ajj  |r(y)  - f(y  + z)J  + o(hM+1) 
and  from  assumptions  (1)  and  (4) 

- o(h"tl) 

Lemma  2 applied  to  Example  #2: 

♦h(v)i  ■ 'x'  [y  * 

♦h(n),  * eX1  [if*  0(h2  )J  . (1  ♦ h)1  (|) 

Therefore 

♦h(v)  + ♦h(n)  " °(h2)»  M " 1 
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Lemma  3. 


From  Lemmas  1 and  2 we  have  Lemma  3. 

Under  assumptions  (1)  through  (4)  of  Lemma  2 and  the  assumption  that  <|>h  Is 
stable  at  Ahy,  It  follows  that  the  solution  nE  to 

♦h(nE)  + s 0 


satisfies 

nE  • a hy  + o(hM+1) 

Indication  of  Proof:  By  Lemma  2 

*h(v) = - *h(ii)  + °(hM+1) 

and  then  applying  Lemma  1 yields  the  desired  result. 

Lemma  3 may  be  generalized  to  theorem  1. 

Theorem  1 ; Assume  that 

(1)  y and  n are  the  unique  solutions  to  F(y)  3 0 and  $h(n)  * 0 where  n Is 

a p-MQ  approximation  to  y with  Mq  ^ p + k and  k ^ 0. 

(2)  Is  a consistent  approximation  to  F of  order  q with  1 < p < q < 2p. 

(3)  and  are  h“^  - Llpschltzlan  at  y. 

(4)  Is  a p-M  approximation  to  F and  <j>h  Is  a Llpschltzlan  p-Nx  approxima- 

tion to  F where  Nx  * min  (M,  MQ  - K,  q - 1 ) . 

Under  the  above  hypotheses.  If  <(>h  Is  stable  at  A^y,  then  the  solution  nE  of 

*h("E)*  ♦h(n)  ■ 0 

satisfies 

nE  ■ V * °(hV‘) 

. 

20 
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SECTION  III 

THE  ERROR  GRADIENT  ESTIMATION  PROCEDURE 

1 . ILLUSTRATION  OF  VARIATION  E ON  ODEs 

In  this  and  the  following  subsection,  we  present  two  methods  or  variations 
for  estimating  the  error  incurred  In  employing  a finite-difference  scheme  to 
obtain  an  approximate  solution  to  an  ODE  initial-value  problem.  For  ease  In 
exposition  we  shall  consider  the  extremely  simple  problem 

= W,  W(0)  = 1 (4) 

along  with  the  difference-scheme 

V-k  » Vk,  V°  = 1 (k  * 0,  1,  2,...)  (5) 

However,  it  will  be  evident  that  the  basic  concept  is  applicable  without  any 
essential  complication  to  more  difficult  problems. 

For  each  value  of  the  parameter  h,  the  mesh-function  Vn  is  filled  In  by 
employing  linear  Interpolation  between  successive  points  (nh,  Vn)  and 
[(n  + 1)  h,  Vn+1];  that  Is,  the  function  U(h,t)  Is  defined  by 

U(h,t)  * Vn(l  - Vn+1  (6) 

where  n * (t/h]  and  t ■ t - nh.  (Recall  that  [x]  denotes  the  largest  integer 
which  does  not  exceed  x.)  In  this  manner  the  function  U(h,t)  Is  defined  for 
all  positive  values  of  h and  nonnegative  values  of  t.  For  fixed  t we  determine 
how  U depends  on  h by  differentiating  equation  6: 

(Since,  for  fixed  t,  the  Index  n undergoes  jump-discontinuities  when  t/h  * 
Integer,  equation  7 must  be  appropriately  Interpreted  for  these  values  of  h as 
a one-sided  derivative,  but  this  does  not  cause  any  difficulty.)  Eliminating 
Vn+1  from  equation  7 with  the  aid  of  equation  5,  one  obtains 
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3U(h,t)  . d\£  / d t \ „n 
3h  dn  l dh  h J v 


In  particular,  as  t/h  -*■  n,  this  equation  assumes  the  form 


aU(h,t)  . dvn  „wn 


The  value  of  V Is  obtained  by  recurrence  from  equation  7,  and  the  value  of 
gjpls  obtained,  also  by  recurrence,  from  the  equation 

dvk+1  - M + hi  dVk  . uk  dV°  n 

dh (1  + h)  dh“+  v • dhT*  0 (1 


which  Is  obtained  by  differentiating  equation  5.  Accepting  the  fact  that 
U(h,t)  approaches  W(t)  as  h 0,  and  assuming  that  Is  integrate  (which 

Is  trivially  true  In  the  present  case  and  is  not  difficult  to  establish  for  any 
first-order  ODE  with  right-hand  side  satisfying  very  mild  hypotheses),  we  obtain 


U(h,t)  - W(t)  » Jh  |^!h,,t) 


In  particular.  If  h = t/n  this  assumes  the  simpler  form 


Vn-W(t).  fhf{|!h,’t>  dh'  (12) 

Jo 

Thus,  one  obtains  the  following  upper  bound  on  the  error  Vn  - W(t): 

|Vn-W(t)|<h  max  |f#h,,t)|  (13) 

0 < h'  < h 3n 

In  the  present  case,  the  explicit  computation  Indicates  that  depends 

monotonlcally  on  h',  so  that  equation  13  becomes  (for  h = t/n): 


|Vn  - W(t)|  < h ||^(h,t)| 


23 


AFWL-TR-78-1 1 9 


In  fact,  a more  accurate  analysis  furnishes  the  improved  estimate 

|Vn  - W(t)|  |f^(h,t)|  (15) 

2.  ILLUSTRATION  OF  VARIATION  H ON  ODEs 

We  consider  the  same  initial-value  problem  and' difference  scheme  as  in  sub- 
section III-l.  The  first  step  Is  to  extend  the  domain  of  definition  of  the 
solution  V of  the  difference  scheme  to  all  nonnegative  values  of  t.  To  do  this 
we  Introduce  the  solution  u(h,t)  to  the  continuum-difference  equation 

u(h»t)  = u(h,t),  u(h,0)  * W(0)  = 1 (16) 

Here,  for  convenience,  we  have  converted  h to  a dimensionless  parameter, 

0 ± h <.  1 , by  Introducing  St  = hAt.  For  0 < t < St  it  is  evidently  necessary 
to  define  u(h,t).  That  Is,  first  we  select  a function  f(h,t)  for  0 <_  t < fit 
(with  f ( h ,0)  * 1)  and  we  define  u(h,t)  for  0 <_  t < <5t  by 

u(h,t)  = f (h,t) 

Then  the  above  difference  scheme  determines  u(h,t)  for  t > St;  clearly  u(h,tn)  = 

Vn  for  h * 1 . If  the  difference  scheme  equation  16  is  convergent,  then  lim 

h+O 

u(h,t)  exists  and  is  readily  seen  to  equal  W(t),  the  solution  of  the 
Initial -value  problem. 

Observe  that  the  solution  to  the  initial -value  problem  is  W(t)  = exp(t); 
to  the  finite-difference  scheme  is  Vn  s (1  + At)n;  and  to  the  continuum 
difference  equation  16  is  u(h,t)  * (1  + 6t)n  f(h,x),  where  n ■ [t/st]  and 
t = t - nst. 

It  Is  evident  that  If  u Is  to  be  continuous  it  is  necessary  to  impose  the 
same  restriction  on  f;  similarly  for  differentiability.  The  most  natural  choice 
of  f is  apparently  f(h,t)  » 1 + t,  0 <_  t < <st,  and  so  we  adopt  the  definition 

u(h,t)  = 1 + t,  0 t < 6t  (17) 
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From  equation  16  we  obtain  difference  equations  for  |f  and  |f,  namely 


|f  (h,t  + fit)  = (1  + fit)  |f  (h.t) 


where 


| f (h,t  + fit) 


R(h,t) 


(1  + fit)  §f  (h,t)  - fit  R(h,t) 


ff  (h,t  + fit)  - u(h,t) 


Furthermore,  from  equation  17  we  obtain 


|f  (h.t)  = 1,  |f  (h.t)  = 0 for  0 < t < fit 


and  let 


It  -fto.t") 


1H  » — (1  tn) 

ah  ah  ' 


Variation  H of  the  error  gradient  estimation  procedure  involves  finding  the 
solution  of  the  following  system  of  difference  equations: 


(1  + fit)  Vn,  V°  = 1 


(1  ♦ at)  |f",  |f°  . i 


(1  + fit)  |f  - fit  Rn 


where 


R"  - §f  ‘ - v" 
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The  error  gradient  estimate  for  variation  H Is 


vn  - w(tn>  * IS 


3.  THEORY  OF  THE  ERROR  GRADIENT  ESTIMATION  PROCEDURE 

dU  3 V 

Since  the  derivatives  ^ and  may  become  discontinuous  at  the  mesh-values 
t = tn,  we  take  only  forward  t derivatives.  If,  similarly,  Is  also  Inter- 
preted as  forward  differentiation,  then 


4r-  V(h,t  * h)  + ")  * lith.t  + h) 


and,  similarly. 


du(h,t  + hat)  . 3u(h,t  + hAt)  . A.  3u(h,t  + hAt) 
TIT  " nr  + At  77 


In  both  the  E and  H variations,  the  error  estimate  Is  based  on  the  identity 


A 

u(hx ,t)  - u(h2,t)  * / f£(h,t) 

J h- 


(under  the  reasonable  assumption  that  u Is  absolutely  continuous  with  respect 
to  h).  If,  as  we  assume. 


w(t)  » 11m  u(h,t) 
h -*■  0 


the  above  Identity  yields 


u(h,t)  - w(t)  ■ j lu, (h'.t > dh, 

Thus,  the  error-estimation  problem  becomes  that  of  estimating  the  Integral 
appearing  In  the  last  equation.  In  the  error  gradient  estimation  procedure, 
the  approximation 


26 


AFWL-TR-78-11 9 


d 


3u  (h\tn) 
ah  * 


dh' 


h 3u.(h,tn) 

n ah 


Is  employed.  Of  course,  the  accuracy  of  this  estimate  depends  critically  on 
the  behavior  of  |jj-  In  the  Interval  [0,h].  In  particular.  If  Is  Increasing 
with  respect  to  h,  the  above  approximation  provides  an  upper  bound  on  the 
error.  When  monotonicity  cannot  be  established,  one  must  settle  for  the 
cruder  upper  bound 


h sup 
0 < h' 


< 


|3u  (h',tn) 

h'ss' 


Returning  to  the  simple  Illustrative  example,  one  sees  that  variation  H 
requires  the  solution  of  a system  of  three  equations,  namely: 


yn+i 

vn+1 

aun+1 

at 


(l  + At)  vn,  v°  =■  l 

n + At>  S".  ff°  - 1 
o + At>  If"  - R"  If0  - 0 


In  contrast  to  this,  variation  E requires  only  the  solution  of  a pair  of 
equations,  namely 


Vn+1  * (1  + At)n,  V°  = 1 


dV 

dh 


n+i 


(1 


+ At)  Ml  + vn  «L. 

dh  v ’ dh 


It  would  therefore  appear  that  variation  H will  Involve  roughly  twice  the  effort 
required  by  variation  E In  order  to  compute  an  error  estimate.  Thus,  at  first 
appearance,  variation  E would  appear  to  be  preferable.  However,  In  variation  E 
the  error  bound  Involves  the  expression  - nVn,  which.  In  the  Illustrative 
example,  works  out  to  -tn  (1  + At)n_1,  which  Is  small  In  comparison  with  the 
separate  terms  jjj—  and  nVn.  That  Is,  the  subtraction  of  two  large  quantities 
which  are  almost  equal  Is  required,  and  this  suggests  that  a noisy  calculation 
may  be  risked.  Thus,  It  Is  difficult  to  say  which  of  the  two  variations  Is, 


H 
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on  the  whole,  preferable  to  the  other;  Indeed,  the  answer  may  depend  critically 
on  the  specific  problem  under  consideration. 

In  the  Illustrative  example  both  variations  lead  to  the  same  estimate, 
namely: 


h f£(Mn)  - - Mil  ♦ 


that  Is,  the  error  estimate  amounts  to 


(1  + At)n  - exp(tn)  « - At(l  + At)n_1tn 


By  employing  the  Identity 


J £u(h,t' ) - w(t')  - w(t')]  dt'  = /7h|^,  ( h ' ,t ' ) dh*  dt 


we  may  make  the  estimate 


N n 

J £u(h,t)  - w(t)]  dt  « h At 


Similarly,  In  the  case  of  PDEs,  the  analogous  estimate 


/•  x . /•  t”  u n 

J J ^u(h;t,x)  - w(t,x)J  dt  dx  ss  h EEff" 
xo  0 j*0  m*0  ^ 


J N 


may  prove  useful,  particularly  when  u may  fall  to  converge  polntwlse  to  w but 
may  converge  to  w In  the  L1  - norm;  similar  estimates  can  be  written  for  the 
case  of  convergence  In  other  Integral  norms. 

4.  ILLUSTRATION  OF  VARIATION  E ON  PDEs 

Paralleling  the  considerations  of  subsections  III-l  and  III-2,  we  shall 
consider  In  the  present  and  the  next  subsection  the  following  problem: 
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§f  + §J-  0,  w(0,x)  - w°(x) 


We  choose  positive  numbers  At  and  ax,  then  define  fit  and  fix,  respectively,  as 
hAt  and  hAx,  where  the  parameter  h satisfies  0 < h £ 1.  Now  consider  the 
difference  scheme 


yn+i  _ yn  vn  yn 


°-  v5  * -°(xj) 


where  x^  * jfix.  Proceeding  by  analogy  with  subsection  III-l,  we  extend  the 
mesh-function  V to  the  entire  half-plane  0<^t<»,  - «<x<«by  defining. 

In  the  rectangle  with  vertices  (n  fit,  jfix),  (nfit,  ( j + 1 ) fix),  ( (n  + 1)  fit,  jfix), 
and  ( (n  + 1)  fit,  (j  +1)  fix),  the  function  U(h,  t,  x)  as  follows  (by  linear 
Interpolation): 


t,  x)  *A  + 8^-+C 


fitfix  + 0 fitfix 


where 


t - nfit,  C * x - j fix,  A * Vj,  B =»  v"+1 


C * vn  - Vn  D * vn+1  - Vn+1 
t VJ+1  Vj,  U Vj+1  Vj 


By  rewriting  equation  19  In  the  form 


vn  + vn 
Vj+i  Vj 


*T  ■ rv3-.  * <’  - 


we  Immediately  observe  that  the  condition  r <_  1 Is  necessary  to  assure  conver- 
gence of  the  difference  scheme  (eq.  19)  as  h -*■  0,  and  It  Is  easy  to  show  that 
this  condition  Is  sufficient  as  well.  In  particular,  the  choice  r = 1 furnishes 
the  trivial  difference  scheme  V*j+1  * vjjj,  reflecting  the  fact  that  the  exact 
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solution  of  equation  18  Is  given  by  w(t,x)  = w°(x  - t);  that  Is,  the  solution 
Is  obtained  by  advancing  the  Initial  data  along  the  lines  Inclined  at  45°  to 
the  x-axIs.  In  order  to  avoid  this  trivial  case,  we  shall  consider  a value  of 
r strictly  less  than  one. 

By  repeated  use  of  equation  21  we  easily  obtain,  for  any  positive  Index  k 


We  now  seek  to  develop  an  expression  for  U(h,  1,  1),  that  Is,  an  approxima- 
tion for  w(l,l).  For  ease  In  exposition  we  choose  r * 1/2,  ax  * 0.1,  At  * 0.05, 
w°(x)  * x;  the  procedure  In  the  gereral  case  will  be  apparent  from  this  simple 
Illustrative  example.  We  choose  a value  of  h slightly  less  than  one,  so  that 
20  fit  < 1 < 21  fit,  10  fix  < 1 < 11  fix.  [That  Is,  the  point  (1,1)  at  which  we 
are  seeking  to  estimate  w(t,x)  lies  In  the  Interior  of  the  rectangle  with 
vertices  at  (20  fit,  10  fix),  (20  fit,  11  fix),  (21  fit,  10  fix),  and  (21  fit,  11  fix).] 
From  equation  22  we  obtain,  upon  making  the  appropriate  replacements, 

WJ  ' 2'k  • (to)  L (?)  (J  ' 1 > <23> 

V ' 1=0 

The  expression  on  the  right  can  be  summed  explicitly,  and  we  obtain 

v?  ' TO  (j  - l)  <2‘> 


Thus,  for  the  values  at  the  vertices  of  the  rectangle  under  consideration,  we 
obtain 


„20  . n t/20  h ..21  h „21  h ..21  h 

V10  °*  V11  To*  V10  • Iff*  V11  “ OT  vn  * W 


From  equation  20  and  the  following  formulas  we  now  obtain 

U(h,  1,  1)  * 0 (26) 

Thus,  In  this  case  the  procedure  of  filling  In  the  rectangles  of  the  mesh  by 
linear  Interpolation  furnishes  the  exact  result,  but  this  Is  due  to  the  simple 
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form  of  the  Initial  data.  If  we  complicate  the  problem  slightly  by  choosing 
w°(x)  * x2,  equations  23  to  26  are  replaced,  respectively,  by 


■s-r*-(vr  5 (!)"•"' 


vj  * iso  (4j2  - 4jk  + k2  + k> 


y20  _ h2  „20  _ 3h2  V21  a llh2  V21  _ 11  h2 

vio  * iff’  vn  ' vio  * isff"»  vn  ‘ m~ 


U(h,  1,  1) 


- 5h2  + 8h  - 2 


Upon  setting  h » 1 In  equation  30  or  In  the  first  part  of  equation  29,  we 
obtain  the  approximation 

w(l,  1 ) ~ 

On  the  other  hand,  upon  differentiating  equation  30  and  setting  h - 1 , we 
obtain 


3U  . 1 

aF  'To 

h»l 


and  so  we  are  led  to  the  error  estimate 


w(l,  1)  - U(l,  1,  1)  « h 


3 U ( h , 1,  1 


In  fact,  since  w(l,l)  ■ 0,  the  true  value  of  the  left  side  Is  - 1/20.  As  In 
subsection  III-1,  we  have  obtained  an  estimate  which  Is  double  the  correct 
value.  We  shall  explain  In  section  IV  how  this  factor  of  2 can  be  removed  In 
the  case  of  the  ODE  example,  and  It  appears  that  a similar  (but  probably  more 
delicate)  argument  could  be  applied  In  the  PDE  case  as  well. 
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5.  ILLUSTRATION  OF  VARIATION  H ON  PDEs 
Consider  the  problem 


ff  + lJ-  °.  w(O.x)  - x°(x) 


(34) 


together  with  the  corresponding  difference  scheme 


n+i  n n n 

vi  - v<  v,-  - v,-  i 

J J + 3 _ q vo 

At  + AX  U*  Vj 


(35) 


where  x^  » jAx.  The  continuum  difference  equation  associated  with  equation  35 
Is  now  written  out: 


u(h»  t + fit,  x)  - u(h;  t,  x)  u(h;  t,  x)  - u(h;  t,  x - fix)  _ p 
5t  fix  ' u 


where  u(h;  0,x)  * w°(x),  u(h;t,x)  * f(h;t,x)  for  0 < t < fit. 

Now  for  h * 1 and  t = tn  (eq.  36)  should  agree  with  equation  35,  and  as 
h -*•  0,  fit  and  fix  should  also  approach  zero.  As  before,  we  take  fit  * hAt.  As 
for  the  relation  between  fix  and  ax,  we  reason  as  follows: 

The  ratio  fix/Ax  must  be  a function  g(h)  such  that  g(l)  and  g(h)  = 0. 
Stability  of  the  difference  scheme  (eq.  35)  requires  that  At/Ax  <_}.  The 
natural  choice  of  g(h)  dictated  by  these  considerations  Is  g(h)  = h.  Then 
equation  36  permits  the  rearrangement 

u(h;  t + fit,  x)  ■ (1  - r)u(h;  t,  x)  + ru(h;  t,  x - fix)  (37) 

where  r * At/Ax.  Equation  37  suggests  how  f should  be  prescribed. 

Let 

u(h;  t,  x)  ■ (1  - r)  w°(x)  + rw°(x  - 

where  0 ^ x < fit  and  5 ■ x/r. 

From  equation  37  we  may  derive  difference  equations 

dU 

we  obtain 
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for  !l  and  Ir-  For 
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It  (h;  t + fit,  *)  • (1  - r)  S-  (h;  t,x)+r^  (h;  t,  x - 6x)  (39) 


d(J 


3U 


3t 


and  from  equation  38  we  obtain 


at 


at 


f£  x)  ■ - (x  - £)  for  0 ^ t < fit 


(40) 


[As  above,  e « r/r  and  0 ^ t < fit.] 
Similarly, 


U (h;  t + fit,  x)  - (1  - r)  (h;  t,  x)  + r (h;  t,  x - fix)  (41) 


au 


ax 


and  again  from  equation  38  we  obtain 


& - «> ■ It  (x)  - ' (l?0(x)  - &°(x  ‘ s) ) *-  o , x < 


fit  (42) 


From  equation  37  we  derive  a difference  equation  for  fjj-,  namely 


ah’ 


§£  (h;  t + fit,  x)  * 0 - r)  (h;  t,  x)  + r|£(h;  t,  x - fix)  - At  R(h;  t,  x) 


(43) 


where 


R(h;  t,  x)  * (h;  t + fit,  x)  + Jjj.  (h;  t,  x - fix) 


au 


Finally,  from  equation  37  we  obtain 


au 

ah 


(h;  t,  x)  ■ 0 for  0 < t < fit 


(44) 


In  variation  H,  the  difference  equations  that  need  to  be  solved  for  the 
error  gradient  estimates  are 

n+i 


v"  ■ 0 • r)  ’J  * v°  ■ “°(xj) 


(45) 
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3U 


n+i 


It  mi  - r) 

j 


£u  + r lu  iy.°  - £w°  / y \ 

at,  '■  ' ' ' atj  r 3tj  " " 3x  (xj) 


dU 


n+i 


. r)  M + r liL  iiL  = aw 

at  '■  r'  a-  T r iw  » 


3t 


— + r — — - — t \ 

3xj  3xj-l’  3xj  " 3X  V j / 


(46) 


(47) 


and 


an  _ (i  . nu  ._n  du  - 

ah  “ ('  " r’  alT  + aTT  ' 4tR.,  tt-  = 0 


n+i 


3un  » 3un 


n 3u 


3h 


j 


3hj  3hj-i  j’  3hj 


(48) 


where 


Rn  aun+1  3un 

J + »J.1 


Note  that  the  amount  of  work  involved  can  be  reduced  by  replacing  equations  46 
and  47  with  the  single  equation 


(49) 


where 


r ■ <’  - r>  ■$  * r RjV  a;-t* 


. (1  . r)  ja  ♦ r 2a 

3tJ  ‘ 


Therefore,  the  difference  equations  to  be  solved  are  reduced  to  equations  45, 
48,  and  49,  and  the  error  gradient  estimate  Is  given  by 


1 j “ xj)  w aKj 


As  a second  Illustration  of  variation  H of  the  error  gradient  estimation 
procedure,  we  consider  the  heat  equation: 


3?  " w(0,x)  ■ w°(x) 


34 


The  difference  equation  may  be  rewritten  In  the  form 


(1  - r)  Vj  + r 


( 


where  r = 2at/ax2.  A familiar  analysis  shows  that  the  tlmestep  restriction 
r £ 1 Is  necessary  and  sufficient  for  stability.  Therefore,  we  let  at  = hat 
and  ax  a h^ax,  and  the  associated  continuum-difference  scheme  Is 

u(h;  t + at,  x)  = (1  - r)u(h;  t,  x)  + r [u(h;  t,  x + 6x)  + u(h;  t,  x - ax)]y 

with  initial  condition 


u ( h ; t,  x)  * (1  - r)  w°(x)  + r [w°(x  + s)  - w°(x  - C)  j/2 


where 


0 < t < 6t  and  c * (2r/r) 


1/2 


Now,  differentiating  with  respect  to  h,  we  obtain 


fr  (h;  t + at,  x ) ■ ( 1 - r)  — r-  (h;  t,  x)  + r 


3U 


3h 


uh 


+ 1^-  (h;  t,  x - ax) 


]/- 


[ft *• 


x + ax) 
at  R (h;  t,  x) 


where 

R(h;  t,  x) 


*7  (h;  t + at,  x)  - 


!“  (h;  t,  x + ax)  - U (h;  t,  x - ax) 
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with  Initial  condition 


(h;  t,  0)  * 0 for  0 < t < 6t 


3h 


Note  that  Rj  * R(1;  tn,  x j ) satisfies  a difference  equation  of  the  form  of 
equation  50.  Therefore,  we  can,  as  before,  replace  the  pair  of  difference 


r)  LI  dU 

equations  for  and  by  a single  difference  equation  for  R. 


Thus,  variation  H leads  to  three  difference  equations: 


n+i 


(1  -r)  vj  *-•(»"«  tVi)/2-,5'"°(xj) 

C - r)^  +0/2‘4t 


and 


>n+i 


(1 


r)  Rj"  ♦ r (nj 


j+1  Rj 


■)/' 


with  the  additional  conditions 


M.'  . (,  . r)  M°  , r tOL°  ) 

nj  8ti  VV  3tj-J 

2JI0 , l!rl(x j ) m l.sai  /„  N 

3tJ  dx*  ’ dXJ  dX  ' J J 


■ ft1  - fe  ♦ f*°  )/• 

3tj  y3xj+i  9 J-yf 


2ax 


In  the  next  example  we  show  why  and  how  f must  be  modified  In  order  to 
obtain  better  estimates  when  the  Initial  data  contain  discontinuities.  Con- 
sider the  problem 


|f  + §f  * °,  w(0,x)  - w°(x) 


i ... - ...  . . 


— 
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The  LKL  (Lax-Keller-leapfrog)  scheme  furnishes  the  following  discrete-difference 
equation: 


vj+1  = 7 (1  ‘ r)  vj+i  + 7 (1  + r)  vj-i 

where  r = At/Ax.  As  before,  the  condition  r < 1 Is  necessary  and  sufficient 
for  stability. 

The  associated  continuum  difference  scheme  Is 

u(h;  t + fit,  x)  = j £u(h;  t,  x + fix)  + u(h;  t,  x - fix ) J 
- £ £u(h;  t,  x + fix)  - u(h;  t,  x - fix )J 


n , 1 


with  Initial  data  u(h;  t,  x)  = f(h;  t,  x)  for  0 x < fit.  If  we  were  to 
proceed  as  in  the  previous  examples,  we  would  define  f by 

f(h;  t,  x)  * j [w°(x  + ?)  + w°(x  - 5)]  - f [w°(x  + c)  - w°(x  - e)] 

for  0 <_  x < fit,  with  % = x/r.  We  refer  to  this  definition  as  the  first  option 
for  f. 


Differentiating  with  respect  to  t we  obtain 


21  (h-  x x)  ■ i-  ( + 5)  _ dw^(x  - eA  1 / dw°(x  + O . dw°l 
3t  in,  x,  x)  -37*  J-2\Zr  +d \T 


(x  - 0 


Next,  letting  x 0 we  obtain 


3f  (h;  0,  x)  a _ dv£(x) 
at  dx 


and,  finally,  letting  h -►  1 and  x -►  Xj  we  obtain 


37 


Similarly,  differentiating  f with  respect  to  u and  letting  t -*•  0,  h -*  1 , and 
x -►  Xj,  we  obtcln 


3u°  . dw°(Xj) 
dx 


Similarly,  differentiating  with  respect  to  h and  carrying  out  the  same  limiting 
operations,  we  obtain 


Note  that  If  = It-  - a 0 for  all  j,  then  f r-  * 0 for  all  Indices  n and  j. 

3hj  anj  3Xj  3hj 
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Now  we  consider  a case  In  which  the  Initial  data  w°  has  a discontinuity. 
In  particular,  let 


w°(x) 


"i  ,f  x * xs 


W If  X > x„ 
r s 


where  w^  and  wr  are  different  constants.  (The  subscript  s suggests  that  a 
shock  Is  being  described.)  If  xs  Is  distinct  from  all  the  numbers  x.,  then 

S U 3 U ft  li  J 

the  Initial  data  for  -rr,  nr,  and  f-  are  zero  at  all  the  x,,  and  so  the 

dn  aun  dX  J 

unacceptable  error  estimate  *>  0 (for  all  n and  j)  would  be  obtained.  To 
remedy  this  situation,  we  Introduce  a modified  definition  of  f,  which  we  term 
the  second  option.  Let 


f(h;  0.  Xj)  - /(xj) 


and 


H(h:  °-  xj)  • ("here  m 


jAX) 


Then  spline-fit  f In  the  Interval  (xj,  x j+1 ) . That  Is,  In  this  Interval  let 
f (h;  0,  x)  be  defined  as  a cubic  polynomial. 


f (h;  0,  x)  = a^x 

where  the  coefficients  aj,  bj, 
the  following  sense: 

f ( h;  0 , Xj  ) 

tx  (h*  °*  xj) 

f(h;  °’  XJ+i) 

ll(h;  °*  xj) 


Xj)3+bj(X-Xj)2+Cj(X-Xj)  + dj 
cj,  dj  are  chosen  to  provide  a smooth  fit.  In 

* w°(xj) 

' ["°(xj«)  ' w°(xj-i)]/(2lx) 

* "0(xj,i) 

■ [w°(xj«)  • “°(xj)]/(2ax) 
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The  coefficients  are  easily  seen  to  be  uniquely  determined  by  the  above  condi- 
tions, and  are  given  by 

*j  ■ (•  “M  *l'j-  Vi  * 7 wj« )/  (lx)3 

bj  ’ (7  "j-i  ' 7 "j  + 7 Vi  ‘ 7 wj«  )/ (4x)2 

CJ  a f£(h*  xj)  ' [w°(xj+.)  ' w°(xj-i)]/  124x1 
dJ  * "°(xj) 

Next,  define  f(h;  t,  x)  for  0 <_  t < 6t  by 

f (h;  t,  x)  = j JV(h;  0,  x + S)  + f(h;  0,  x - £)] 

- § [f(h;  0,  x + c)  - f(h;  0,  x - 5)] 

where  £ * x/r.  Taking  the  t-derlvatlve  of  f and  letting  x -►  0,  one  obtains 


|£{h;  0,  x)  * - |£(h;  0,  x) 


Now  let  h -*■  1 and  x -*■  Xj  In  order  to  obtain  the  discrete  Initial  conditions  for 

3u  3u°  . 

at*  it  j 


Similarly,  differentiating  f with  respect  to  x and  letting  x -*•  0,  h -*■  1 , and 
x -*■  Xj  furnishes 

- * <=. 

9XJ  J 

Note  that  f Is  actually  Independent  of  h,  therefore 

du°  Q 

V 
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! 


! 


1 


wm 


Observe  that  the  second  option  for  f,  In  contrast  to  the  first,  does  not 
lead  to  Identically  vanishing  Initial  values  for  |£  and  ~ when  the  Initial 
condition  Is  given  by  a step  function  (actually  containing  one  or  more  discon- 
tinuities). Therefore,  the  second  option  Is  to  be  preferred  when  discontinuous 
data  are  prescribed. 

An  additional  advantage  of  the  second  option  Is  that  the  program  user  Is 
not  required  (as  In  option  1)  to  Input  for  In  option  2 the  program  computes 
[w°(xj+i)  - w°(xj_1 )]/(2ax)  from  the  given  w°(xj).  Thus,  with  option  2 the 
modified  program  (l.e.,  modified  to  furnish  error  estimates)  does  not  require 
any  Information  beyond  that  required  by  the  unmodified  program. 

The  next  example  Illustrates  the  use  of  the  procedure  In  the  case  of  a non- 
linear PDE.  Consider  the  problem 

|f  * f£(w)  - 0,  „<0,x)  ■ «°(x) 

where  F Is  assumed  to  possess  suitable  smoothness  properties.  The  LKL 
scheme  furnishes  the  discrete-difference  equation 


vn+1 

Vj 


where  r - At/Ax  and  Fj  * F(vj)’  The  e(luat1on  °f  first  variation  Is 


.n+i 


*.  3J-*  . r /An  6n  _ An  n \ 
2 i rj+i  j+i  Aj-i  6j-i/ 


where  Aj  * F'(vj)*  This  can  be  rearranged  to 


+ r 


U 


where  rj  * rAj.  A stability  analysis  shows  that  the  condition  |r!j|  <_  1 Is 
necessary  for  stability.  Therefore,  we  take  fit  * hAt  and  fix  » hAx. 


I 
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The  associated  continuum-difference  scheme  Is 

u(h;  t + st,  x)  * j[u(h;  t,  x + 6x)  + u(h;  t,  x - x)J 

- £ [F(h;  t,  x + 6x)  - F(h;  t,  x - fix)] 

where 

F(h;  t,  x)  * F[u(h;  t,  x)] 

The  initial  condition  Is  given  by 

u(h;  t,  x)  = f (h;  x,  x)  for  0 <_  x < fit 
We  choose  the  second  option  for  f: 

f(h:  0-  Xj)*  “°(Xj) 

fj‘hi  0’  Xj)  ’ [w°(v)  ' w°(xJ-i)]/(2Sx> 

and  then  we  spline-fit  f,  for  x between  x^  and  xj+1 , by  a cubic  polynomial, 
exactly  as  In  the  previous  example.  Next  we  define 

f ( h ; t,  x)  for  0 < x < fit  by  f (h;  x,  x)  - j [f(h;0,  x + £)  + f(h;0,  x - o] 

“ 7 [^(h;  x + 5)  - F°(h;  x - c)] 

where  C * x/r  and  F°(h;  x)  * F[f(h;  0,  x)].  Now,  forming  ~ and  letting  x -*•  0, 
we  obtain 

If h;  °* x)  * ■ §7° (h;  x)  8 ■ A[f<h;  °. *)]  ff  o,  x) 
where  A ■ F' . Next,  let  h ■*  1 and  x Xj*,  one  obtains 

42 


Similarly,  by  differentiating  f with  respect  to  x and  letting  t «►  0,  h -*■  1 , 
x -*>  Xj , we  obtain 


Also,  differentiating  f with  respect  to  h and  letting  t -*■  0,  h -*■  1 , x -*>  x j 
we  obtain 


Observe  that  the  foregoing  procedure  extends  to  the  case  when  u,  F,  and  f 
are  vectors  (In  which  case  A becomes  a matrix). 
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SECTION  IV 

NUMERICAL  EXPERIMENTS 

1 . INTRODUCTION 

The  method  of  asymptotic  error  estimation  and  the  error-gradient  method 
were  tested  experimentally  on  four  ODEs,  In  this  section  we  present  the 
numerical  results  In  both  tabular  and  graphical  form, 

2,  THE  ODEs  AND  OAEs 

In  analogy  with  the  standardized  abbreviation  ODE  for  ordinary  differential 
equation  we  shall  employ  the  abbreviation  OAE  for  ordinary  difference  equation. 

The  ODEs  considered  here  are  of  the  form  dw/dt  * F(w).  Four  specific  cases 
are  discussed,  namely:  F(w)  * w,  F(w)  = - w2,  F(w)  « 0.9  w2,  and  F(w)  = - w3. 

In  each  of  these  cases  the  exact  solution  can  be  written  explicitly,  and  this 
enables  us  to  determine  the  true  values  of  the  errors  Incurred  by  the  difference- 
scheme  and  thus  to  compare  the  accuracy  of  the  several  methods  employed  to 
estimate  these  errors.  In  this  subparagraph  IV-2,  the  exact  solutions  and  the 
output  of  the  difference-schemes  are  presented  for  each  of  the  four  equations; 
in  subparagraph  IV-3,  the  estimates  furnished  by  the  asymptotic  error  estimation 
procedure  are  worked  out;  while  In  subparagraph  IV-4,  the  estimates  obtained  by 
the  error  gradient  procedure  are  obtained. 

In  each  case  we  take  for  the  OaE  the  Euler  scheme: 

Vn+1  * Vn  + At  F(Vn) 

The  values  of  Vn  are  termed  the  exact-OAE  solution.  Of  course,  the  actual 
output  of  the  computer  Is  only  an  approximation  to  the  values  of  Vn.  Certain 
OaEs  will  be  generated  whose  solutions  serve  as  estimates  of  the  error  Incurred 
by  using  Vn  as  an  approximation  to  the  exact  solution  of  the  ODE.  These  are 
the  error  estimates,  whose  values  are  approximated  by  the  computer  output. 

Of  course,  the  real  challenge  Is  to  obtain  reliable  error  estimates  for  ODE 
and  PDE  problems  whose  exact  solutions  are  not  available;  nevertheless,  use  of  a 
proposed  error-estimation  procedure  on  some  problems  with  known  solutions  is 
frequently  employed  to  furnish  an  Indication  of  how  well  the  procedure  will  work 
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In  more  challenging  cases.  Furthermore,  the  analysis  of  simple  cases  often  turns 
out  to  be  valuable  In  debugging  the  programs  that  are  written  to  obtain  the 
desired  output. 

Case  1 

- w,  w(0)  - 1 


The  solution  of  the  ODE  problem  Is  given  by  w * et,  while  the  solution  of 
the  OaE  problem  ^namely,  Vn+1  * (1  + At)Vn,  V°  * ij  Is  given  by  Vn  ■ (1  + At)". 

Remark;  It  is  of  Interest  to  note  that  there  exists  In  this  case  a 
difference  scheme  which  yields  the  exact  solution,  namely  Vn+1  ■ eAt  Vn.  When 
such  an  exact  difference  scheme  Is  known  It  can  be  useful  In  analyzing  both  the 
error  and  the  error  estimate.  This  Is  Illustrated  In  paragraph  3,  In  the  dis- 
cussion of  the  error  estimate  furnished  by  the  asymptotic  error  estimation  pro- 
cedure for  Case  1 . 

Case  2 

- - w2,  w(0)  * 1 

The  exact  solution  of  the  ODE  problem  Is  given  by  w ■ 1/(1  + t).  The  OaE 
problem  Is 


vn+1  « (l  - A tv" ) v",  V°  - 1 

While  the  computational  scheme  Is  quite  trivial,  the  exact  solution  of  the  OaE 
problem  cannot  be  expressed  In  simple  form.  However,  one  can  obtain  a partial 
representation  by  the  following  device. 

Replacing  At  by  -k  and  then  replacing  kVn  by  Wn,  we  convert  the  above  OaE 
Into  the  parameter- free  form 

Wn+1  - (l  + w")  w"  (51) 

with  Initial  condition  W°  - k.  Thus,  W1  - k + k2,  W2  - (k  + k2)  + (k  + k2)2  - 
k + 2k2  + 2k3  + k4,  W3  * (k  + 2k2  + 2k3  + k4)  + (k  + 2k2  + 2k3  + k4)2  • k + 3k2 
+ 6k3  + 9k4  + 10k5  + 8k6  + 4k7  + k8,  etc.  By  Induction,  It  Is  clear  that  Wn  Is 
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a polynomial  of  degree  2n  In  k,  with  positive  Integer  coefficients  (except  for 
the  vanishing  constant  term).  Symbolically,  we  can  write 


Wn  ■ a^,  *()_!»••••  3j,  where  N * 2n 


Then  the  recurrence  relation  (eq,  51)  furnishes  a recurrence  relation  for  the 
coefficients  a”,  namely, 

an+1  - an  + T a"  a" 
m m **  j m-j 

J 

where  It  Is  understood  that  a^J  * 0 when  m < 1 and  when  m >_  N.  This  Is  what  we 
mean  by  a partial  representation;  It  Is  useful  for  debugging  purposes. 

Remark:  It  Is  Interesting  to  note  that  the  exact  solution  to  the  OaE 
problem 


Is  given  by 


0n+I  - Dn(l  ♦ kDnw).  0°  ■ 1 


Dn  * (1-nk)’1 


Upon  replacing  k by  -At  we  obtain  Dn  * (1  + nAt)-1,  which  agrees  exactly  with 
the  solution  of  the  ODE  problem  under  discussion  here.  Thus,  as  In  Case  1,  a 
slight  modification  of  the  Euler  scheme  turns  out  to  furnish  the  exact  solution. 

Case  3 


- 0.9  w2,  w(0)  - 1 


The  exact  solution  of  the  ODE  problem  Is  given  by  w ■ 1/(1  - 0.9  t)  for 
0 £ t < 1 , 

Vn+1  - (1  +0.9  AtVn)Vn,  V°  « 1 
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As  In  Case  2,  no  convenient  representation  of  Vn  Is  available;  of  course,  the 
entire  discussion  of  the  preceding  case  Is  applicable  If  certain  minor  changes 
are  carried  out. 

Case  4 


g^  = - w3,  w(0)  * 1 

The  exact  solution  of  the  ODE  problem  Is  given  by  w * (1  + 2t)_Js.  The  OaE 
problem  Is 

Vn+1  = (l  - At  (V")2)  V".  V°  * 1 

A partial  representation  of  the  solution  of  the  OAE  problem  can  be  obtained  by 
following.  In  a rather  obvious  manner,  the  above  discussion  of  Case  2. 

3.  ASYMPTOTIC  ERROR  ESTIMATES 

In  this  section  a brief  review  Is  given  of  the  asymptotic  error  estimation 
procedure  for  ODEs,  and  then  numerical  results  for  each  of  the  four  specific 
cases  listed  In  paragraph  2 are  presented.  Finally,  a typical  FORTRAN  program 
for  performing  the  required  computations  Is  given, 

a.  Review  of  the  Procedure  for  ODEs 

Given  the  ODE  g-j:  * F(w),  let  <t>  and  <(>E  be  two  finite-difference 
approximations  of  the  operator  g£  - F,  <f>E  being  of  higher  order  than  $.  For 
example, 

n vn+1  »n 

♦ (V)"  * -£~- -V-  - F(Vn) 

and 


are  first-order  and  second-order  approximations,  respectively.  Let  V be  the 
solution  to  «p( V)  * 0 and  let  V£  be  the  solution  to  ♦(v^)  + <frE(V)  ■ 0,  each 
satisfying  the  prescribed  Initial  condition.  (The  standing  assumption  Is  made 
that  the  solutions  V and  V£  exist  and  are  unique.)  Then  V - VE  Is  defined  to 
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be  the  asymptotic  error  estimate,  That  Is,  V - V£  is  taken  as  an  approximation 
to  be  actual  error  V - w Incurred  by  using  the  lower-order  scheme  <t> ( V ) * 0, 


For  a given  choice  of  $ there  is  wide  freedom  In  choosing  the  higher- 
order  approximation  to  - F.  In  order  to  estimate  what  effect  the  latitude 
In  choosing  $£  may  actually  have  In  practice,  we  have  considered  two  specific 
choices  In  our  numerical  experiments: 


1 


Option  A: 


i 

Option  B: 


4A)(v)n 

4B)(V)n 


While  the  lower-order  approximation  <t>  was  chosen  as  above,  namely 

♦(«*  - - f(v") 

The  solution  of  the  lower-order  scheme  <f> ( V)  = 0 is  simply  denoted  as  V,  while 
the  solutions  of  the  two  higher-order  schemes  are  denoted  and  V^;  that 
Is, 


♦(V)  * 0 

* (4A>)  + 4A)(V)  - 0 

* (4B)) + 4B)(v)  = 0 

From  the  theory  of  asymptotic  error  approximations  (ref,  1),  the  difference 
V-VE*V-w  + o(at2),  [Therefore,  * o(At2),]  What  our 

numerical  experiments  will  provide  Is  some  Indication  of  the  magnitude  of  the 
constants  Implied  In  the  o(At2)  terms  for  each  of  the  options  A and  B, 

In  order  to  give  a heuristic  overview  of  the  asymptotic  error  estimation 
theory,  we  Illustrate  how  It  works  In  Case  1,  where  all  the  calculations  can  be 
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. 


worked  out  quite  explicitly.  [Recall  that  F(w)  * w,  w(0)  ■ 1.]  As  pointed  out 
previously,  In  this  case  there  exists  an  exact  difference  scheme,  namely  Vn+1  * 
eAt  Vn;  thus,  If 

„n+i  Atun 

mv)"  , *-.%•  v 

the  difference  operator  a Is,  In  an  obvious  sense,  an  exact  discretization  of 
D * - (•);  that  Is,  the  operator  - Identity.  The  existence  of  an 

exact  difference  scheme,  while  exceptional,  is  helpful  In  Illustrating  just 
what  Is  occurring  In  the  calculations. 

Now,  let  be  an  oUt*)  approximation  to  D;  then 


*k(V)N  * Mv)n  - kJ  Atk 


where  Kk  * 0(1),  a bounded  quantity.  In  particular,  let 
fk(At)  = 1 + fy  + + •••  + (k  = 

and  let 


„n+i  vn 

,k(V)"  . 1—^1  . fk(4t)Vn 


Since 


♦k(v") 


/ eAt  - 1 - 4tfk(4t)  \ 

■ "<*>  + \ !t— 5 /v 


It  follows  that  Is,  In  fact,  an  o(Atk)  approximation  to  D,  with 


"Kk  Atk  * ( 


eAt  - 1 - Atfk(At) 
At 


)vn 


or 


Kk  * kk(V)  * °0)Vn  + o(At) 
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From  the  above  formulas  we  Immediately  confirm  the  following  observation: 

Lemma  A:  For  k * 1,  2,  3,  •••  the  quantity  k£(V)  satisfies  a Llpschltz 
condition  modulo  o( At ) ; that  Is,  there  exists  a positive  constant  £ such  that 


iK^Vj)  - Kj(V2)|  < £|vj  - v!J| 


Two  additional  observations  appear  worthy  of  note.  The  proofs,  which 
are  rather  straightforward,  are  omitted. 

Lemma  B:  If  <|>k  Is  an  o(Atk)  approximation  to  D,  then  the  actual  error 
Is  o(Atk).  That  Is,  let  ^(V)  * 0 and  D(w)  - 0;  then  Vn  - w(tn)  * o(Atk), 

Lenina  C.  Let  k and  i be  positive  Integers,  k < and  let  <f>k  and  <j>  be 
o(Atk)  and  o(At*)  approximations,  respectively,  to  D.  Let  w be  the  solution  to 
D(w)  = 0,  let  V be  the  solution  to  $k(V)  = 0,  and  let  be  the  solution  to 
4>k(VE)  + 4>£(  V)  = 0.  Then 


Vn  = w(tn)  + o(Atk) 

v£  = w(tn)  + o(At2k)  + C^t*1) 

Vn  - v£  = Vn  - w(tn)  + o(At2k)  + o(At£) 


( 52-a ) 


(52-b) 


(52-c) 


If  In  Lemma  C we  take  l ■ 2k,  then  the  three  conclusions  assume  the 


»"  - «(t")  * o(atk) 

VJ  - w(tn)  ♦ o(At2k) 
vj  - Vn  - w(tn)  + o(At2k) 


(53-a) 


(53-b) 


(53-c) 


Thus,  In  Case  1,  It  becomes  very  clear  why  Vn  - v|?  constitutes  an  estimate  of 
Vn  - w(tn)  of  order  2 k 

Now,  as  Indicated  above,  we  present  some  numerical  results  for  each  of 
the  four  problems  listed  In  subparagraph  IV-2. 
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Case  1 

In  this  case,  options  A and  B coincide: 


♦ E(V) 


n _ Vn+1  - yn  Vn+1  + Vn 


and  therefore 


v£+1  - 0 + At)  v2  + (At)2  Vn/2 


By  Iteration  It  Is  easily  shown  that  this  equation  (together  with  the  initial 
condition  = 1 ) Implies 


v£  = Vn  + (1  + At)0"1  . tnAt/2 


Therefore,  the  exact-OAE  error  estimate  produced  by  the  asymptotic  error  estima- 
tion procedure  In  the  present  case  (for  either  option)  Is  given  by 


Vn  - vj?  * -(1  + At)"'1  • tnAt/2 


n-i  j. 


while  the  actual  error  Is  Vn  - enAt,  or  (1  + At)n  - enAt,  which  Is  equal  to 


and  so 


-enAt  . tn  At/2  + o( At2 ) 


vn  - Vp  * Vn  - enAt  + o(At2) 


Choosing  At  * 0.1  and  n * 10  we  obtain:  w(l)  * e « 2.718, 

V10  * (1.1 )10  « 2.594,  and  V^°  « 2.712.  The  actual  error  Is  therefore 


V10  - w(l)  * - 0.124 
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while  the  asymptotic  error  estimate  Is 

V10  - V£°  a - 0.118 

Upon  reducing  At  to  0.05  and  correspondingly  Increasing  n to  20,  we  obtain 
w(l)  * e a 2.718,  V20  a 2.653,  and  V|°  a 2.716.  Thus,  the  actual  error  is  now 

V20  - w(l)  a - 0.065 


while  the  asymptotic  error  estimate  becomes 


V20  - V20  a - 0.063 


Case  2 

In  this  case,  option  A reduces  to 


4A)(V)n  = (Vn)4  (At)2/4  - (V")3At 


while  option  B assumes  the  form 


4B)(V)n  = (Vn)4  (At)2/2  - (Vn)3At 


As  In  Case  1,  we  first  choose  At  = 0.1  and  n * 10,  then  At  = 0.05 
and  n * 20.  With  the  first  choice  we  obtain  w(l ) * 0.5,  V10  a 0.48171, 
^EA  * 0*49976,  and  a 0.49943,  Thus,  the  actual  error  is 

V10  - w(l ) % - 0.01829 

while  option  A furnishes 


V10  • a - 0,01805 
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and  option  B furnishes 


V10  - V|g  * - 0.01772 

With  the  second  choice  of  At  and  n we  obtain  w(l)  * 0.5,  V20  « 0.49110, 
V2°  « 0.49994,  and  V|g  « 0.49986,  and  so 

V20  - w(l)  w - 0.00890 
V20  - Vgjj  « - 0.00884 
V20  - V|[j  « - 0.00876 


Case  3 

With  At  = 0.1  and  n * 10  we  obtain  w(l)  * 10,  V10  « 4,46663, 
VEA  * 9<08783>  V|g  « 9.20051,  Thus,  the  actual  error  Is 

V10  - w(l)  « - 5.53337 

while  option  A furnishes  the  error  estimate 

Vl°  “ vea  * ’ 4-60040 


and  option  B yields 


vl°  - veb  * " 4-73388 

When  At  Is  reduced  to  0.05  and  n Is  Increased  to  20,  we  obtain 
w(l)  - 10,  V2°  w 5.72982,  V|j  « 13.80567,  V2°  * 13.92144.  Thus,  the  actual 
error,  option  A estimate,  and  option  B estimate  are,  respectively, 
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V20  - w(l)  a - 4.27018 
V20  - V|0  » - 8.07585 
V20  - V2°  a - 8.19162 


The  rather  violent  behavior  of  the  output  Is  perhaps  to  be  expected  In 
view  of  the  rapid  growth  of  the  exact  solution.  However,  the  fact  that  the 
error  estimates  obtained  with  At  » 0.05  are  worse  than  those  obtained  with 
At  * 0,1  Is  somewhat  unexpected. 


Case  4 


With  At  = 0.1  and  n » 10,  we  obtain  w(l)  = 3"1*  » 0.57735,  V10  « 0.56044, 
^EA  * 0-57679,  « 0.57628.  Thus,  the  actual  error  Is 


V10  - w(l)  a - 0.01691 


while  option  A furnishes  the  error  estimate 


vl°  - V^°  a - 0.01635 


and  option  B yields 


V10  - V*°  a - 0.01584 


When  At  Is  reduced  to  0,05  and  n Is  Increased  to  20,  we  obtain 
w(l)  » 0.57735,  V20  * 0.56917,  V|j|  * 0.57722,  V|jj  as  0.57709. 

Thus,  the  actual  error,  option  A estimate,  and  option  B estimate  are, 
respectively. 


V20  - w(l)  a - 0.00818 


V20  - V|»  a - 0.00805 


V20  - V|j  « - 0.00792 
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b.  FORTRAN  Program  for  Asymptotic  Error  Computations 

For  option  A,  with  At  = 0.1,  the  required  computations  were  carried  out 
with  the  very  simple  FORTRAN  program  which  appears  In  figure  1.  Reduction  of 
At  to  0.05  required  only  a change  In  line  11,  while  option  B required  only  very 
simple  changes  In  lines  19  through  22. 


I PROGRA F ASYMAP1C INPUT, OUTPUT) 



Vi=i.O 

*2=1.0 

5 V 3=1 . 0 

— — — V4=1.0  — 

W10LD=1. 0 

weoto=iv  o 

W30L0=1. 0 

-to— W4OL0-1.  a — — 

M = 0 • 1 

23  PRINT  2 0,  T . VI.  V2.  t/3,  V4*  H10L0  « W20L0  .W30L  0,  W40L0 

20  FORMAT  <F5.2»8F15.  10) 

W1  W1  etf H-H*  Kioto 

15  W2NEW=  l*20L  0«H*W20L0*'*2 

W3N£  W=W3OLO«-0.  9*H*W30L0**2 

W4N£W=W40Lt)-H*W40 L0**3 

T=T*M- 

Vl=Vl«-H*Vl-WlNtW*H10L0*HM  W10L0*WlNEW)/2. 0 

-20 *2=*2-H*  V2*“*2-W2  N€  W2  Ot  O-H* 1 f M2  Ofc  (H  H2  N£ W > /2 . 0 > * * 2 

V 3= V3*H* V3*#2- W3N£ W+W33LD+  0. 9*H* ( ( W30L0+N3NEH)  /2.  0)  **Z 
V4«V4-H*V4**3-M,N£W4-W40L0-HM  < M40LG+H4NEM) /2 . 0 ) **3 
W10L0= N1NEH 
*»20L0=  W2NE  W 
25  W3OL0=  H3NEM 

- H I,  A«  H/'LI  

WHULu*  H^rtc  — * 

IFU.LT. 1.0)60  TO  23 

STOP 

ENG 


Figure  1.  FORTRAN  Program  for  Asymptotic  Error  Computations 
4.  ERROR  GRADIENT  ESTIMATES 

We  begin  this  section  with  a very  brief  review  of  the  error  gradient 
procedure,  which  Is  the  subject  of  section  III.  For  simplicity  In  exposition, 
the  ODE  problem  will  be  assumed  to  be  of  the  form  ^ ■ F(w),  w(0)  ■ 1 , and  the 
finite-difference  technique  to  be  employed  will  be  taken  as  the  Euler  method. 
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vn+1  « vn  + hF(Vn) , V°  - 1 


(54) 


The  discrete  function  V,  with  values  V°,  V1,  V2,***  Is  then  filled  In,  In  the 
Initial  Interval  [0,h],  more  or  less  arbitrarily,  by  a function  U(h,t)  which 
satisfies  U(h,0)  - V°  - 1 and  U(h,t)  * V1.  Then  the  definition  of  U(h,t)  for 
t > h Is  accomplished  by  using  the  obvious  extension  of  equation  54,  namely 


U(h,t  + h)  - U(h,t)  + hF  (u(h,t)) 


For  any  given  value  of  t (for  simplicity,  t Is  chosen  equal  to  1 In  the  lllus- 
tratlve  computations  which  follow),  ^ Is  determined  as  a function  of  h. 

[if  U(h,t)  lacks  sufficient  smoothness,  must  be  suitable  Interpreted.] 

Then,  for  any  given  (positive)  value  H of  h,  one  obtains 


U(H,1)  - 11m  U(h,l ) * /*H  U(h,1)  dh 
h*o  3 


Under  very  mild  hypotheses  the  limit  appearing  In  this  equation  exists  and 
equals  w(l),  the  value  of  the  actual  solution  at  t * 1 . Thus,  the  error  Is 
given  by 

H 


w(l)  - U(H,1)  * - f |^(h»1}  dh 

ft 


While  a rigorous  justification  has  not  been  obtained,  the  estimate 


"f;11 


(55) 


has  furnished  remarkably  good  results  In  the  four  cases  which  are  considered  In 
this  report. 

In  each  of  the  four  problems  given  In  paragraph  2,  three  definitions  of 
U(t,h)  were  employed  In  the  Interval  0 <_  t <_  h.  First,  U(t,h)  was  defined  to 
be  linear  In  this  Interval;  second,  U(t,h)  was  defined  to  be  quadratic  In  the 
Initial  Interval,  the  three  available  coefficients  being  so  chosen  that 
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U(h,o)  * V°  ■ 1 , U(h,h)  * V1,  and  remains  differentiable  when  t * h 

(and  therefore  at  all  succeeding  nodes);  and,  third,  U(t,h)  was  defined  to  be 
cubic  In  the  Initial  Interval,  the  four  available  coefficients  being  chosen  so 
that,  In  addition  to  the  three  constraints  Imposed  In  the  quadratic  case, 
U(h,t)  possesses  a continuous  second  derivative  at  the  nodes. 

Thus,  12  computations  were  carried  out  altogether;  In  each  case  U(h,l)  was 
determined  for  h * 10  J'V  k - 1 » ^ L n i.  100,  1 <_  k <_  10,  so  that  h assumed 
the  values  ygjj,  ygy,  yjjj-,  •••  . The  quantities  the  error 

estimates  - h and  the  corrected  estimates  U( h ,1 ) - h were 

also  computed;  this  final  quantity  was  also  plotted.  The  12  plots  (figures  3 
through  14)  thus  obtained  are  Included  at  the  end  of  this  report.* 

We  list  here  a few  of  the  numerical  results  (tables  1 through  4);  however, 
the  plots  give  a very  clear  picture  of  the  manner  In  which  the  corrected  esti- 
mates approach  the  actual  solution  as  h -*■  0. 

c.  FORTRAN  Program  for  Error  Gradient  Computations 

The  FORTRAN  program  for  Case  2 (w1  * -w2,  w(0)  = 1)  with  quadratic 
Interpolation  appears  In  figure  2;  this  program  provides  for  plotting  the 
corrected  estimate  of  the  solution,  namely,  U(h,l)  - h h » 1 ) . father 
straightforward  modifications  of  this  program  were  used  for  the  other  11 
computations. 


I 
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♦The  captions  of  the  plots  are  taken  from  the  titles  of  the  computer  programs; 
CASE1L,  CASE1Q,  ....  CASE4C.  The  numeral  Indicates  which  of  the  four  differ- 
ential equations  Is  being  solved,  while  the  final  letter  refers  to  linear, 
quadratic,  or  cubic  Interpolation, 
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TABLE  1 


U (h,l)  - h||  - w(l)  f 

CASE  1 [w(l ) = e 

as  2.718282]- 

h 

CASE  1L* 

CASE  IQ 

CASE  1C 

1 

T5 

.111255 

-.012256 

-.010388 

l 

10.5 

-.013551 

-.008648 

-.009326 

1 

IT 

(-.114083 

-.r’0292 

-.008722 

j .102934 

l 

Iff 

(-.064984 

-.003351 

-.002838 

j .061363 

1 

2TTT5 

-.003897 

-.002462 

-.002707 

1 

2T 

(-.062019 

-.003053 

-.002585 

j .058721 

1 

(-.013603 

-.000147 

-.000125 

j .013444 

i 

995 

-.000179 

-.000112 

-.000123 

1 

TOT 

(-.013468 
j .013312 

-.000144 

K II 

-.000122 

♦In  the  case  of  linear  Interpolation  |^-  Is  discontinuous  when  h Is  the 
reciprocal  of  an  Integer;  the  left-hand  and  right-hand  limits  are 
therefore  listed. 


, 
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TABLE  2 


U(h,l)-h  |^(h,1) 

- wtl)  , CASE  2 

[w(l)  - 0.5] 

h 

CASE  2L 

CASE  2Q 

CASE  2C 

1 

IS 

.026860 

.002606 

.001040 

1 

io.  5 

.000918 

.000220 

.000932 

l 

n 

j -.024940 
{ .024246 

.002113 

.000846 

i 

19 

013133 
) .012937 

.000588 

.000239 

i 

293" 

.000226 

.000061 

.000227 

1 

2T 

j-. 012477 
{ .012301 

.000531 

.000216 

1 

97 

j-. 002550 
] .002542 

.000022 

.000009 

1 

997 


.000009 


.000003 


,000009 


1 

nnr 


.002524 

.002517 


.000022 


,000009 
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TABLE  3 


0(h>1)  - hfh  - 

w(l)  , CASE  3 [w(l) 

- io] 

h 

CASE  3L 

CASE  3Q 

CASE  3C 

1 

10 

-2.674869 

-3.869900 

-3.825690 

1 

105 

-3.744252 

-3.699798 

-3.720758 

1 

IT 

j -4 . 679396 
| -2.477286 

-3.661614 

-3.621111 

1 

20 

) -3 .387269 
j -1.320952 

-2.398601 

-2.376689 

1 

20.5 

-2.341970 

-2.317923 

-2.328550 

1 

21 

j -3.284377 
j -1.236982 

-2.302749 

-2.282015 

1 

9? 

(-.896927 
| .168615 

-.368956 

-.366563 

1 

99.5 

-.365858 

-.362621 

-.363807 

1 

100 

(-.888060 
] .170631 

-.363436 

-.361083 
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TABLE  4 

Ll(h .1 ) - w(l),  CASE  4 [w(l)  - 3’1*  » 0.577350] 


h 

CASE  4L 

CASE  4Q 

CASE  4C 

1 

10 

.031889 

.004662 

.001207 

1 

10.5 

.001264 

-.000510 

.001067 

k-  ! 

11 

1 

030592 

| -028680 

.003723 

.000970 

1 i 

20 

i 

-.015641 

.015101 

.000976 

.000262 

1 

20.5 

.000306 

-.000092 

.000248 

1 

21 

i-. 014836 

1 .014349 

.000879 

.000237 

1 | 
77 

i-. 002961 

1 .002941 

.000035 

.000010 

l 

99.5 

.000-012 

-.000003 

.000010 

loo  j 

1 -.002931 

1 .002911 

.000034 

.000010 
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SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 

Admittedly,  the  computations  on  the  error-gradient  method  reported  here  are 

a 

much  more  extensive  than  those  on  the  ^symptotlc  method.  In  Cases  1,  2,  and  4, 
the  results  obtained  with  the  asymptotic  method  are  better  than  those  obtained 
by  using  the  error-gradient  method  with  cubic  Interpolation,  which  In  turn  Is 
superior  to  the  error-gradient  method  with  either  linear  or  quadratic  Inter- 
polation. However,  In  Case  3,  reduction  of  the  step-size  from  h * 1/10  to 
h * 1/20  results  In  a distinct  worsening  of  the  asymptotic  results,  while  all 
three  Interpolation  methods  furnish  Improved  values.  All  results  In  Case  3 are 
disappointing,  but  this  Is  to  be  expected  In  view  of  the  rapid  growth  of  the 
solution. 

A notable  phenomenon  Is  the  bracketing  effect  obtained  In  all  four  cases 
when  linear  interpolation  Is  employed  (see  figures  3,  6,  9,  and  12  for  Cases  1L, 
2L,  3L,  4L,  respectively).  Except  In  Case  3,  where  Irregular  behavior  Is  to  be 
expected,  the  bracketing  sets  in  from  the  very  beginning  (h  * 1/10),  and  the 
upper  and  lower  envelopes  are  essentially  linear  with  Intercept  at  the  exact 
solution.  In  Case  3 oscillation  appears  Immediately,  but  bracketing  does  not 
appear  until  h has  diminished  to  0.02,  approximately.  Hwould  certainly  be 
helpful  If  one  could  prove  that,  when  linear  Interpolation  Is  employed,  the 
error-gradient  method  furnishes  values  which  bracket  the  true  solution  for 
sufficiently  small  values  of  h. 

It  Is  also  of  Interest  to  note  that  when  linear  interpolation  Is  employed 
and  the  left-  and  right-hand  limits  are  averaged  when  h Is  the  reciprocal  of  an 
Integer,  a very  marked  Improvement  Is  obtained.  Thus,  In  Case  1,  the  average 
error  when  h * 1/20  Is  « - 0.0018,  while  V|°-  w(l)  « - 0.002  (see  paragraph  3); 
similarly.  In  Case  2 the  average  error  when  h * 1/20  Is  »-  0.00020,  while  the 
errors  - w(l)  and  - w(l)  are  approximately  - 0.00024  and  - 0.00057, 
respectively.  In  Case  4,  the  three  corresponding  numbers  are  - 0.00027, 

- 0.00056,  and  - 0.00107.  Closely  related  to  these  observations  Is  the  fact 
that  when  h Is  the  reciprocal  of  a half-integer  (e.g.,  1/20.5)  the  error  Is 
noticeably  dlmlshed;  e.g.,  h ■ 1/20,5  gives  an  error  much  smaller  than  the 
errors  obtained  either  with  h ■ 1/20  or  h ■ 1/21, 


! 


I 


63 


AFWL-TR-78-1 1 9 


It  Is  by  no  means  clear  that  a satisfactory  theory  can  be  worked  out,  but  It 
certainly  appears  that  further  numerical  experimentation  and  theoretical  Investi- 
gations are  justified. 
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